home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / RANDMZJ.C < prev    next >
C/C++ Source or Header  |  1996-03-26  |  5KB  |  175 lines
  1. /* ============ */
  2. /* randmzj.c    */
  3. /* ============ */
  4. #include <defcodes.h>
  5.  
  6. /* ------------------- */
  7. /* FUNCTION PROTOTYPES */
  8. /* ------------------- */
  9. # undef F
  10. # if defined(__STDC__) || defined(__PROTO__)
  11. #    define  F( P )  P
  12. # else
  13. #    define  F( P )  ()
  14. # endif
  15.  
  16. /* INDENT OFF */
  17. static    long    DEK_rand F((void));
  18. extern    int    IRand_MZ F((void));
  19. extern    long    LRand_MZ F((void));
  20. extern    double    Rand_MZ F((void));
  21. extern    void    Rand_MZ_Init F((unsigned long));
  22. extern    void    VRand_MZ F((double *, int));
  23.  
  24. # undef F
  25. /* INDENT ON */
  26.  
  27. #define TWOM24    1.0/TWOP24        /* 2^(-24) */
  28. #define TWOP24    16777216.0        /* 2^24    */
  29. #define P24MASK 16777215L        /* 2^24-1  */
  30. /* INDENT OFF */
  31. static    double    Seeds[24] =
  32. {
  33.      4138572.*TWOM24, 16412963.*TWOM24,  5060164.*TWOM24,  3954602.*TWOM24,
  34.     11488643.*TWOM24, 11099381.*TWOM24, 10909270.*TWOM24,   564186.*TWOM24,
  35.     14261907.*TWOM24, 10716489.*TWOM24, 10042238.*TWOM24, 14709425.*TWOM24,
  36.      7010471.*TWOM24, 15189851.*TWOM24,  1616330.*TWOM24, 10844364.*TWOM24,
  37.      8026029.*TWOM24, 13538854.*TWOM24,   983044.*TWOM24,  3024518.*TWOM24,
  38.     14475343.*TWOM24,  5142119.*TWOM24,  4187065.*TWOM24,   764925.*TWOM24
  39. };
  40. /* INDENT ON */
  41.  
  42. static double Carry = 0;
  43. static int i24 = 24, j24 = 10;
  44. /* ==================================================================== */
  45. /* Rand_MZ - Returns uniform variate by method of Marsaglia-Zaman-James    */
  46. /* ==================================================================== */
  47. /* ------------------------------------------------------------ */
  48. /* The source document for this generator is:            */
  49. /*                                */
  50. /*    F. James, A review of pseudorandom number generators,    */
  51. /*    Computer Physics Communications 60(1990) 329-344.    */
  52. /* ------------------------------------------------------------ */
  53. double
  54. Rand_MZ(void)
  55. {
  56.     double  NextRand;
  57.  
  58.     NextRand = Seeds[i24 - 1] - Seeds[j24 - 1] - Carry;
  59.  
  60.     if (NextRand >= 0)
  61.     {
  62.     Carry = 0;
  63.     }
  64.     else
  65.     {
  66.     ++NextRand;
  67.     Carry = TWOM24;
  68.     }
  69.  
  70.     Seeds[i24 - 1] = NextRand;
  71.  
  72.     --i24;
  73.  
  74.     if (i24 <= 0)
  75.     {
  76.     i24 = 24;
  77.     }
  78.  
  79.     --j24;
  80.  
  81.     if (j24 <= 0)
  82.     {
  83.     j24 = 24;
  84.     }
  85.  
  86.     return NextRand;
  87. }
  88. #define FULL_SIZE ((unsigned)RAND_MAX + 1U)
  89.  
  90. /* --------------------------------------------- */
  91. /* IRand_MZ - Returns Standard C Uniform Variate */
  92. /* --------------------------------------------- */
  93. int
  94. IRand_MZ(void)
  95. {
  96.     return((unsigned) (Rand_MZ() * (double) FULL_SIZE) & RAND_MAX);
  97. }
  98. /* ----------------------------------------------- */
  99. /* LRand_MZ - Returns Uniform Variate of Type Long */
  100. /* ----------------------------------------------- */
  101. long
  102. LRand_MZ(void)
  103. {
  104.     return((long) (Rand_MZ() * TWOP24));
  105. }
  106. /* ------------------------------------------------------- */
  107. /* VRand_MZ - Returns Vector of Uniform Variates in [0, 1) */
  108. /* ------------------------------------------------------- */
  109. void
  110. VRand_MZ(double *RandVec, int VecLen)
  111. {
  112.     int     k;
  113.  
  114.     for (k = 0; k < VecLen; ++k)
  115.     {
  116.     RandVec[k] = Rand_MZ();
  117.     }
  118. }
  119. static unsigned long RandSeed;
  120. /* ------------------------------------------- */
  121. /* DEK_rand - From Knuth vol 2, p 102, line 26 */
  122. /* ------------------------------------------- */
  123. static
  124. long    (DEK_rand) (void)
  125. {                    /* Knuth Random Numbers */
  126.     /* Knuth Generator Returns a 24-Bit Integer */
  127.     RandSeed = RandSeed * 1664525 + 1;
  128.     return((unsigned long) ((RandSeed >> 8) & P24MASK));
  129. }
  130. /* ------------------------------------------------------------------ */
  131. /* Rand_MZ_Init - Initializes Marsaglia-Zaman Random Number Generator */
  132. /* ------------------------------------------------------------------ */
  133. void
  134. Rand_MZ_Init(unsigned long NewSeed)
  135. {
  136.     int     k;
  137.  
  138.     RandSeed = NewSeed;
  139.  
  140.     /* --------------------- */
  141.     /* Warm up the generator */
  142.     /* --------------------- */
  143.     for (k = 1; k <= 100; ++k)
  144.     {
  145.     DEK_rand();
  146.     }
  147.     for (k = 1; k <= (sizeof(Seeds) / sizeof(double)); ++k)
  148.     {
  149.     long    NextRand = DEK_rand();
  150.  
  151.     Seeds[k - 1] = (double) NextRand *TWOM24;
  152.     P(printf("Seeds[%2d] = %.15e\n", k + 1, Seeds[k]));
  153.     }
  154. }
  155. # if defined(TEST_RAND_MZ)
  156. void
  157. main()
  158. {
  159.     int     k;
  160.  
  161.     Rand_MZ_Init(341055109L);
  162.  
  163.     for (k = 1; k <= 50; ++k)
  164.     {
  165.     double    NextRand = Rand_MZ();
  166.     printf("NextLongRand # %2d = %10.1f\n", k, NextRand * TWOP24);
  167.     }
  168.     for (k = 1; k <= 50; ++k)
  169.     {
  170.     int    NextRand = IRand_MZ();
  171.     printf("NextIntRand  # %2d = %5d\n", k, NextRand);
  172.     }
  173. }
  174. # endif
  175.